This section presents a mathematical overview of surface parameterization, B-Spline surface representation and least squares based B-Spline surface approximation techniques.
A. Surface Parameterization
Recently a great number of 3D models are available and can be easily processed to triangle meshes. However, triangle mesh representation of free-form surfaces, although flexible, suffers from its high dimensional complexity (three coordinates per vertex) in a recognition point of view, and therefore the hardness of manipulation. Surface parameterization, as a technology to construct a mapping between the meshes and some parameter domain, therefore to reduce the 3D problems to 2D ones, is well studied and discussed.
Surface parameterization calculates specific kinds of mapping from surface itself to some parameter domain, e.g., the 2D plane. Surveys like provided in M. S. Floater and K. Hormann, “Surface parameterization: a tutorial and survey,” in In Advances in Multiresolution for Geometric Modelling. Springer, 2005, pp. 157-186 and in A. Sheffer, E. Praun, and K. Rose, “Mesh parameterization methods and their applications,” in Foundations and Trends in Computer Graphics and Vision. Now Publishers, 2006, pp. 105-171, give a comprehensive introduction to this topic.
In an application of reverse engineering, generally speaking, due to the spirit of tensor product B-spline, or NURBS, the parameter domain is a rectangle, or a polygon whose edges are all parallel to the axes. According to Equation (2), the spline surface can be acquired by first evaluating the spline curve in each of the x, y dimensions respectively, and then using tensor product structure to build the regular gird (see FIG. 9 for an example). It is highly desired that the iso-parametric lines are locally orthogonal to reduce the conversion error. To achieve this, Harmonic Mapping is employed to compute such a mapping.
Any mapping (u(x; y); v(x; y)) which satisfies the following two Laplace equations is called a harmonic mapping.
A conformal mapping f satisfies the Cauchy-Riemann equations, which, with z=x+iy and w=u+iv, are
                    ∂        u                    ∂        x              =                  ∂        v                    ∂        y              ,                    ∂        u                    ∂        y              =                            ∂          v                          ∂          x                    .      By differentiating one of these equations with respect to x and the other with respect to y one obtains the two Laplace equations:Δu=0,Δv=0,where
  Δ  =                    ∂        2                    ∂                  x          2                      ⁢                  +                  ∂        2                    ∂                  y          2                    is the Laplace operator.
In the continuous domain, the parameterization problem can be formulated as follows: given a genus O surface S⊂R3 with one boundary, compute a map f:f:S→R2,f(∂S)=Γwhere Γr is a closed 2D curve.
There are infinitely many mapping from a surface to a 2D disk with fixed boundary. According to M. O. Carmo, Differential Geometry of Curves and Surfaces. Prentice Hall, 1976, the mapping from a surface embedded in R3 to R2 will cause distortions, either angle distortion, or area distortion, or both, unless the surface is developable. Here some method is sought to minimize the Dirichlet functional, or the string energy:
      E    ⁡          (      f      )        =            1      2        ⁢                  ∫        S            ⁢                                                              ∇              f                                            2                .            This mapping is usually regarded as a harmonic map, which has a low angle distortion and is easy to compute.
The energy optimization is equivalent to solving the following PDE:Δf=0,  (1)where Δ is the Laplacian operator.
B. B-Spline Surfaces
B-Splines are piecewise polynomial functions. A tensor product B-Spline surface is given by:
                              S          ⁡                      (                          u              ,              v                        )                          =                              ∑                          i              =              0                        m                    ⁢                                    ∑                              j                =                0                            n                        ⁢                                          Q                ij                            ⁢                                                N                                      i                    ,                                          d                      0                                        ,                                          t                      0                                                                      ⁡                                  (                  u                  )                                            ⁢                                                N                                      j                    ,                                          d                      1                                        ,                                          t                      1                                                                      ⁡                                  (                  v                  )                                                                                        (        2        )            where m and n are the number of control points in the u and v parametric directions. Ni,dl,tl, l=1, 2 . . . are B-Spline functions of degree dl with knot vector tl. QijεR3 are control points of the surface. See E. Cohen, R. Riesenfeld, and G. Elber, Geometric modeling with splines: an introduction. AK Peters Ltd, 2001 for further details about B-Splines and tensor product B-Spline surfaces, which is sufficiently known to one of ordinary skill.
C. Least Squares Surface Approximation and Variants
Given a regular grid of 3D points, a tensor product B-Spline surface can be fit to approximate the data grid such that the distance of the approximated surface and the data points is minimized. Least squares based surface approximation is widely used for this task as for instance described in E. Cohen, R. Riesenfeld, and G. Elber, Geometric modeling with splines: an introduction. AK Peters Ltd, 2001, and [7] L. Piegl and W. Tiller, The NURBS book. Springer Verlag, 1997.
The objective function to be minimized is given by
                                          F            lsq                    ⁡                      (            Q            )                          =                              ∑                          k              =              0                        r                    ⁢                                                                                    S                  ⁡                                      (                                                                  u                        k                                            ,                                              v                        k                                                              )                                                  -                                  P                  k                                                                    2                                              (        3        )            where Pk are the sample data points, (uk,vk) are the parameter values of Pk and Qi,j are the unknown variables. The degree is appropriately set; for e.g., cubic in both parametric directions.
Flsq is minimized by setting its partial derivatives with respect to Qi,j to zero and solving the resulting linear system for the unknown control points. The knot vectors are computed such that the every knot interval has at least one sample. This ensures that the resulting linear least squares matrix is positive definite as is known to one of ordinary skill and is also described in L. Piegl and W. Tiller, The NURBS book. Springer Verlag, 1997 and C. De Boor, A practical guide to splines. Springer Verlag, 2001.
1) Weighted Least Squares:
The least squares technique may be enhanced by setting weights to samples to increase or decrease their influence for surface fitting. The least squares objective function may be modified as follows.
                                          F            wlsq                    ⁡                      (            Q            )                          =                              ∑                          k              =              0                        r                    ⁢                                    w              k              2                        ⁢                                                                                                S                    ⁡                                          (                                                                        u                          k                                                ,                                                  v                          k                                                                    )                                                        ⁢                  _                  ⁢                                                                          ⁢                                      P                    k                                                                              2                                                          (        4        )            where wk is the assigned weight for Pk.
2) Smoothing Least Squares Splines:
The least squares techniques tend to introduce extraneous unwanted wiggles. These wiggles are small scale high curvature features. In order to obtain smooth surfaces, the least squares objective function may be modified to bound the magnitudes of second derivatives of the surface that in turn bounds the magnitude of the curvatures. The modified objective function is given by:
                                          F            smooth                    ⁡                      (            Q            )                          =                                            ∑                              k                =                0                            r                        ⁢                                          w                k                2                            ⁢                                                                                                            S                      ⁡                                              (                                                                              u                            k                                                    ,                                                      v                            k                                                                          )                                                              -                                          P                      k                                                                                        2                                              +                      λ            ⁢                                          ∑                                  k                  =                  0                                r                            ⁢                              (                                                                                                                                                    S                          uu                                                ⁡                                                  (                                                                                    u                              k                                                        ,                                                          v                              k                                                                                )                                                                                                            2                                    +                                                                                                                                    S                          uv                                                ⁡                                                  (                                                                                    u                              k                                                        ,                                                          v                              k                                                                                )                                                                                                            2                                    +                                                                                                                                    S                          vv                                                ⁡                                                  (                                                                                    u                              k                                                        ,                                                          v                              k                                                                                )                                                                                                            2                                                  )                                                                        (        5        )            where Suu, Suv, Svv are the second derivatives of the surface S(u,v). It is to be noted that the resulting system obtained after setting partial derivatives of Fsmooth to zero is still linear since the smoothing term is also a quadratic function of the control points. The smoothing parameter (λ) influences the quality of the approximation and must be chosen judiciously to obtain satisfactory results.
3) Constrained Least Squares Splines:
It is possible that some control points are fixed by external constraints. In this case the unknowns are the remaining control points and a surface approximation can be obtained by formulating a constrained least squares objective function. First rewrite the tensor product surface as follows.
      S    ⁡          (              u        ,        v            )        =            ∑              t        =        0            M        ⁢                  Q        l            ⁢                        N          l                ⁡                  (                      u            ,            v                    )                    M=mn,Ql=Qi(l)j(l) Nl(u,v)=Ni(l)(u)Nj(l)(v)i(l)=[l/n],j(l)=l mod n Qa, a⊂l=0 . . . M are fixed control points. Let L=l−a. The constrained least squares objective function is given by equation (6):
                                          F            cons                    ⁡                      (            Q            )                          =                              ∑                          k              =              0                        r                    ⁢                                                                                                        ∑                                          l                      =                                              L                        0                                                                                    L                                                                      L                                                                                                      ⁢                                                            Q                                              l                        ^                                                              ⁢                                                                  N                                                  l                          ^                                                                    ⁡                                              (                                                                              u                            k                                                    ,                                                      v                            k                                                                          )                                                                                            -                                  (                                                            P                      k                                        -                                                                  ∑                                                                              a                            ^                                                    =                                                      a                            0                                                                                                    a                                                                                  a                                                                                                                          ⁢                                                                        Q                                                      a                            ^                                                                          ⁢                                                                              N                                                          a                              ^                                                                                ⁡                                                      (                                                                                          u                                k                                                            ,                                                              v                                k                                                                                      )                                                                                                                                )                                                                    2                                              (        6        )            
Workflow
In this section each step of the executed pipeline, from images to meshes, spline surfaces and finally CAD models in details is provided.
A. Image Acquisition and Mesh Generation
The high resolution CT image is first acquired. And some segmentation method is applied to extract the region of the bone from the background in each images. After that, a mesh generation method like Marching Cubes which is known to one of ordinary skill and is for instance described in W. Lorensen and H. Cline. “Marching cubes: A high resolution 3d surface construction algorithm,” in Proceedings of the 14th annual conference on Computer graphics and interactive techniques. ACM. 1987. p. 169 gives the triangle mesh of the model.
B. Mesh Processing and Preparation
After the mesh of the entire model is obtained, it is further processed to get a decent one as the input of parameterization. First, the region of interest is highlighted and a cut a generated along the boundary between the region and other parts of the model. Then the mesh of highlighted region is extracted, with a special care of the quality of the cutting boundary.
Due to the consideration of numerical stability, smoothing or remeshing is performed (e.g. the algorithm as described in M. Botsch, M. Pauly, L. Kobbelt, P. Alliez, B. Levy, S. Bischoff, and C. Rossi, “Geometric modeling based on polygonal meshes,” in SIGGRAPH '07: ACM SIGGRAPH 2007 courses. New York, N.Y., USA: ACM, 2007) if necessary. The resulting mesh should not contain very thin triangles, or floating triangles (ones two of whose edges are boundary edges).
C. Mesh Parameterization
Earlier herein, harmonic map was briefly introduced and why it is fitted to the present purpose. Next the implementation details are explained.
A triangle mesh M consists of a vertex set V, an edge set E, and a face F. Each vertex v in V is given its position in the Euclidean space R3, i.e., v={xv, yv, zv}. The connection of vertices, or the topological information is encoded in the edge set and face set. Each edge e in E is connecting two vertices in V, i.e. eij=[vi, vj]. Similarly, three consecutive loop edges form a face in F, i.e., fijk=[vi, vj, vk]. In the rest of this document, sometimes p is used instead of v to represent the vertices to avoid the confusion of similar symbols.
Using a finite element method, U. Pinkall and K. Polthier, “Computing discrete minimal surfaces and their conjugates,” Experimental Mathematics, vol. 2, pp. 15-36, 1993 was able to find a solution to the discrete counterpart of Equation by defining the discrete Laplacian operator on the triangle mesh. It is known as the cotangent formula. The discrete Laplacian off at each vertex vεV in mesh M is:
      Δ    ⁢                  ⁢          f      ⁡              (        v        )              =            ∑                                    v            i                    ∈          V                                      [                          v              ,                              v                i                                      ]                    ∈          E                      ⁢                  (                              cot            ⁢                                                  ⁢                          α                              v                ,                                  v                  i                                                              +                      cot            ⁢                                                  ⁢                          β                              v                ,                                  v                  i                                                                    )            .      αv,vi and βv,vi and are defined as shown in FIG. 10. When [v, vi] is a boundary edge, only one angle is considered.
Setting the Laplacian to 0 yields the solution to the harmonic map. Usually, such a process involves solving two linear systems, one for the x coordinate, and the other for the y coordinate. Program 1 shows in pseudo-code how to calculate the harmonic map from a disk-topology surface to a unit square. Note that the corner points, or the boundary condition need to be chosen carefully to reduce the distortion.
Program 1: Calculate the harmonic map from a disk to a unit squareRequire: Triangle Mesh M with genus 0 and one boundaryEnsure: Δƒ(v) = 0 for v ∉ ∂M , and ∂ƒ is a unit squareSegment ∂M into four consecutive segments Γi (0≦i≦3)Mark all vertices unfixedfor each vertex p in Γ0 doSet u(p) = 0, and mark p as fixedend forfor each vertex p in Γ2 doSet u(p) = 1, and mark p as fixedend forfor each vertex p in M doif p is not fixed thenAdd Δu(p) = 0 to the linear system    end ifend forSolve linear system     to get u(p) for each vertex pMark all vertices unfixedfor each vertex p in Γ1 doSet v(p) = 0, and mark p as fixedend forfor each vertex p in Γ3 doSet v(p) = 1, and mark p as fixedend forfor each vertex p in M doif p is not fixed thenAdd Δu(p) = 0 to the linear system    end ifend forSolve linear system    to get v(p) for each vertex pƒ(p) = (u(p), v(p))
The parameterization from a disk to other planar shapes is essentially the same. The only difference is to set different boundary conditions for each linear system.
D. Sample Grid Creation
The mesh is then resampled based on the parameterization result. For the square case, see Program 2 in pseudo-code for details. The template case is basically the same.
Sometimes the distortion along the boundary is much higher than other region. This can be solved by extracting a little more areas in the preprocessing step, then drop the boundary of the regular grid.
Program 2: Uniform sampling of the mesh to get a regular grid:Require: Parameterized Triangle Mesh M with u(p) and v(p) for each vertex p, and sampling resolution RES for i ← 0 to RES do  for j ← 0 to RES do   In the parameter domain, find triangle Δ which con-         tains    ⁢                  ⁢    planar    ⁢                  ⁢    point    ⁢                  ⁢          p      ij        =      (                  i        RES            ,              j        RES              )     Calculate the Bari-centric coordinate of pij with re-   spect to Δ   Apply the coordinate interpolation in the 3D domain   to get a sampling point Sij  end for end for Connect Sij to form a regular grid
E. B-Spline Surface Fitting
This section presents two methods for tensor product B-Spline surface fitting using least squares based techniques presented in an earlier section. First, an iterative method is presented for fitting a single B-Spline surface to a rectangular grid of samples with associated parameter values. Second, an algorithm is presented for fitting a multi-patch surface taking cross-boundary smoothness into consideration. Bi-cubic surfaces have been used for all results, and appear satisfactory. The iterative methods are variations of the approach presented in V. Weiss, L. Andor, G. Renner, and T. Varady, “Advanced surface fitting techniques,” Computer Aided Geometric Design, vol. 19, no. I, pp. 19-42, 2002.
1) Single Patch Fitting: Given a rectangular grid of samples, an accurate and sufficiently smooth B-Spline surface is required. Least squares fitting without the smoothing term gives accurate surfaces but contain unwanted wiggles due to the nature of the fitting. The wiggles are not part of the original data set and cause additional problems when creating solid models. At the same time, a very large smoothing parameter will not give accurate surfaces. Therefore smoothing parameter values must be judiciously selected. FIG. 11 shows a surface fit with plain least squares fitting. This surface contains many unwanted small high curvature features. FIG. 13 shows a surface fit using a large smoothing parameter value. This surface is not a sufficiently accurate fit to the sample grid. FIG. 12 shows a surface obtained in accordance with a method provided in accordance with an aspect of the present invention. The surface is smooth as well as accurate. Comparison of the effects of smoothing parameter on least squares fitting is demonstrated in FIGS. 11-13. The control net size (53×53) and knot vectors for all three surfaces are the same. FIG. 11 No smoothing, λ=0; FIG. 12 Optimal smoothing, λ=3.4*10−7; FIG. 13 Over smoothing, λ=10−4.
FIG. 14 shows the variation of the least squares residual error with increasing smoothing parameter value. Each curve is labeled with a number that indicates the number of B-Spline control points (same in U and V). FIG. 15 shows the variation of the maximum fitting error with increasing smoothing parameter value. Each curve indicates variation for a given number of control points as labeled in the legend. The maximum fitting error is computed as the maximum of the distance of a sample to its closest point on the approximated surface. The two graphs are closely co-related. A more detailed analysis is provided below.
First, when smoothing is not used, the maximum error and the least squares residual error progressively decrease as the number of control points increases. This is expected since there are more degrees of freedom to approximate the shape of the data grid. However, as previously noted, such surfaces may have undesirable wiggles and the number of wiggles increases with control net size. Second, for any given control net size, the error growth with increasing smoothing parameter value can be distinguished into three regions. Initially, the error does not increase much. Then there is a point after which the error grows exponentially (linear in graph since it is shown on logarithmic scale). Finally, after a given value the accuracy of the surface is no longer preserved and the surfaces assume minimal curvature. At every iteration, the method estimates the value of the smoothing parameter a little to the right of end of the first segment by exponential interpolation. At this value, the surface will be smoother while at the same time will retain accuracy. Therefore, the goal of the iterative method is to determine the control net size as well as the desired smoothing parameter such that the maximum error bound is satisfied.
As an aspect of the present invention an iterative method is provided using least squares fitting with smoothing to automatically compute optimal number of control points, location of control points as well as smoothing parameter values. The surface fitting method is designed to reduce user assistance as much as possible. The user needs to specify two parameter values—one for the maximum allowed error and one for the amount of smoothness required. At each iteration, the knots are determined so that there is at least one sample within every knot interval. This is done so that the least squared system is positive definite. The method proceeds in the two phases:
Phase I: Determine minimum number of control points (Cmin) such that surface is within maximum error bound using least squares approximation without smoothing. This is implemented using a binary search on the number of control points in the range degree+1 through min(mr,mc).
Phase II: Determine minimum number of control points such that surface is within maximum error bound and satisfies smoothness criterion using least squares approximation with smoothing. The smoothness condition for a given number of control points is deemed satisfied if the total error grows only by a user specified amount. This phase is implemented using a binary search on the number of control points in the range Cmin through min(mr,mc). For a given number of control points, the desired smoothing weight is iteratively estimated using exponential interpolation. The method stops when the maximum error bound is within the tolerance.
Typical values for the maximum allowed error are 0.01 or 0.001. From experiments it was determined that allowing a growth of the least squares residual by one or two orders of magnitude seems to give sufficiently smooth surfaces. FIG. 16 illustrates a knee implant region mesh extracted from a femur. FIG. 17 illustrates a smooth and accurate surface fit related to FIG. 16 using a method as provided herein as an aspect of the present invention. It has a 53 by 53 control net and with λ=3.4×10−7.
FIG. 17 shows the optimal surface computed using the above algorithm for a distal femur data set shown in FIG. 16. An alternative iterative algorithm that successively adds degrees of freedom instead of performing binary search is presented elsewhere herein.
For the data set shown in FIG. 16, the samples are relatively sparse near the front due to the distortion in parameterization. Surface quality may be additionally improved by performing weighted least squares.
FIGS. 18 and 19 compare surfaces fit with and without weights, wherein FIG. 18 shows a surface without least square weights (64 by 64 control net, λ=8.5×10−9) and FIG. 19 shows a surface with least square weights (53 by 53 control net, λ=3.4×10−7). In this example, samples that are closer to the medial and lateral condyle ends were assigned weight value 10.0 since the parametric distortion in these regions results in sparser samples in the grid. FIG. 19 shows that one can get better quality surfaces by setting appropriate weights.
2) Multi-patch Fitting: The distal femur surface has a highly curved outer boundary with concave regions. For such surfaces, a rectangular parameter domain may result in large relative distortions as shown herein. This results in irregular sample grids that will not allow accurate B-Spline surface approximation in regions where sampling is sparse. In general, a polygonal template domain that more closely conforms to the shape of the outer boundary will result in less parametric distortion and therefore more regular sample grids. For the distal femur bone, a template illustrated in FIG. 20 has been chosen. This template domain gives better parameterization results than a rectangular domain as shown in FIG. 27. For surface fitting, the samples are divided into three patches with anatomical correspondences; D1 on left and D3 on right corresponding to condyles and D2 in center corresponding to the patellar surface. Surfaces are fit to the three patches so that C(2) continuity is obtained across shared patch boundaries. The chosen template also has particular advantages for distal femur surfaces since it enables creation of implants for partial knee replacements. This would involve fitting surfaces for only one of the condyles with or without the patellar surface.
The iterative method provided herein as an aspect of the present invention may be used to fit surfaces separately to patches D1 and D2. Also provided as an aspect of the present invention is another method that successively adds degrees of freedom where required. The left and right patches have simpler geometry as compared to the surface as a whole. As a side note, the successive approach when applied to the single patch surface took longer to complete since many more iterations were required.
Recall that in the iterative method provided herein, binary search is performed on the control net size in both phase I and II. The knot vectors at each step were determined based on the control net size and the distribution of samples. This algorithm starts from the minimum control net size (4×4 for a bicubic surface) and successively adds degrees of freedom (knots) in the regions of maximum error. If the same locations are requested multiple times, knots are added in neighboring regions to avoid placement of multiple knots that reduce B-Spline continuity. Knots are placed such that every knot interval has at least one sample.
Surfaces fitted separately to all three patches will not guarantee continuity and smoothness at shared patch boundaries. Therefore after fitting surfaces to patches D1 and D3 separately, a constrained surface is fitted to the central patch DD. C(0), C(1) and C(2) continuity conditions at the shared boundaries of D2 constrain the placement of the first three and last three control vertex columns of D2.
The number of control vertex rows and knots in the vertical direction for D2 are dictated by patches D1 and D3. The knots for the shared boundaries of D2 are refined to match D1 and D3. Additional knots and control vertices are added to satisfy a maximum error tolerance using an iterative algorithm that uses constrained least squares fitting in each step.
FIG. 21 shows an optimal multi-patch surface fit to template based samples. Patch D1 has 9×9 control points, patch D2 has 7×6 control points and patch D3 has 6×7 control points. The multi-patch surface has an order of magnitude fewer control points than the single patch surface that has 532 control points. FIG. 22 shows a map of the Gaussian curvatures of the multi-patch surface. Smooth variation of the curvatures across shared boundaries confirms C(2) smoothness of the surface. FIG. 23 shows a magnified view of the map in the region around the central patch for further visual clarity.
F. CAD Model Creation
In order to convert a B-Spline surface to a thick solid model, an offset surface is first created along the surface normal direction. Isocurves corresponding to the boundaries of the original surface as well as the offset surface are extracted.
Ruled surfaces are created between corresponding pairs of boundary curves. The ruled surfaces are then stitched with the original and offset surfaces along respective common boundary curves. Features such as pins required to assemble the implant onto a patient's bone may be modeled separately and added to the implant solid model using Boolean operations. Other assembly features include screw threads.
FIGS. 24 and 25 illustrate solid models of implants created from B-spline surfaces, with FIG. 24 being a single patch surface with pins of which a pin 2401 is identified and FIG. 25 a multi-patch surface.
FIG. 24 shows solid models created from single and multipatch surfaces obtained from the iterative fitting algorithm. The solid models are then exported to standard CAD formats such as IGES or STEP that can then be sent to a CAM system for manufacturing. Solid model creation, assembly feature addition and CAD format export was performed using OpenCascade as described in S. Open Cascade, “Open cascade technology ver 6.3,” 2010, URLwww.opencascade.org.
Results from Other Data Sets
Next some results of the parameterization are shown. The methods described were tested on an Intel Pentium D 2.8 GHz PC with 2 GB RAM. They were implemented in C++ and compiled by Microsoft Visual Studio 2008. Matlab R2008a was used as the linear system solver. In all tests where the maximum number of vertices of the models is approximately 20,000, the running time is less than 2 seconds. Later in this section the pros and cons of the different parameterization methods for the application of reverse engineering will also be analyzed.
In FIG. 26 the parameterization result of a hearing aid 2601 to a unit square 2602 is shown. On top right in 2602 a (u, v) coordinate for each vertex of the original mesh is assigned, and thus flatten the model to a planar square. On bottom left that coordinate is used as the texture coordinate, a check-board as the texture, and employ texture mapping. From the figure, one can see that the intersection angle of iso-parametric lines is nearly 90 degree, which demonstrates the low stretch, or angle distortion of the mapping. However, as stated earlier, if the mapping is angle-preserving, then usually there should be area distortion, or metric distortion. On the bottom right of FIG. 26 a further mapping illustrates this.
The above method has also been tested on a knee bone model in a similar manner as described above. Due to the apparent concavity of its boundary, it was also tried to map it to an irregular planar. FIG. 27 illustrates the results of a knee bone model. The first row in 2701 shows the original triangulation. The second and third rows show the planar domain, the texture mapping, and the map of the area distortion of two different parameter domains respectively, similar to what has been done in FIG. 26. FIG. 28 demonstrates the low distortion of the multipatch parameterization. FIG. 29 shows a polygonal mesh of a distal femur bone obtained after processing a CT scan image. FIG. 30 shows a single surface patch fit to the mesh after parameterization and grid sampling using the techniques presented herein in accordance with an aspect of the present invention. FIG. 31 shows a template based multi-patch surface fit to the mesh. The outer boundary of the multi-patch surface seems to have sharper corners. This is due to the fact that the mesh parameterization generated degeneracies at the corners due to the poor quality of the mesh triangulation in these regions. Therefore a small area of the mesh near the boundary was chopped and not used for surface fitting. Such degeneracies may be avoided by remeshing the data to obtain better quality triangulations and reparameterizing the new meshes.
In summary, as an aspect of the present invention a complete workflow for creating CAD models of patient-specific customized implants from scanned medical images has been provided. One implementation demonstrates the medical images. One implementation demonstrates the feasibility of creating a unified system to perform all tasks in the proposed pipeline viz., image processing, discrete surface creation and processing, smooth surface creation and solid model creation and processing. A unified system as provided herein enables faster turnaround times for creating customized implants of better quality, thereby improving the market potential of a product based on it.
Anatomical objects typically have highly curved complicated surface geometry. In particular, techniques for mesh parameterization, sampling, single and multi-patch surface fitting, and solid model creation were implemented. The techniques have been demonstrated for the application of creating distal femur implants for knee replacements. The distal femur surface is highly curved, has a highly curved outer surface boundary with concave regions resulting in two finger like extensions (condyles) and a central connecting region (patellar surface). The results indicate that the techniques provided herein as one or more aspects of the present invention are well suited for the problem.
The template based parameterization and successive multipatch surface fitting in particular seems to have potential for fitting other anatomical objects with complicated geometry.
In one embodiment of the present invention, appropriate template shapes are selected for each type of anatomical object and the parameterization and surface fitting techniques are extended for the chosen templates.
A system illustrated in FIG. 32 and as provided herein as an aspect of the present invention is enabled to receiving, processing and generating data. The system is provided with data that can be stored on a memory 1801. Data may be obtained from a data source, for instance via an input 1806. Such data may be medical image data or any other image data. The data may also be CAD data such as a CAD model of a medical implant. The system has a processor 1803. The processor 1803 is provided or programmed with an instruction set or program executing the methods of the present invention that is stored on a memory 1802 and is provided to the processor 1803, which executes the instructions of 1802 to process the data from 1801. Data, such as image data or CAD data or CAM data such as a CNC file or any other signal resulting from the processor can be outputted on an output device 1804, which may be a display to display data or a data storage device. The output device may also be a storage device that stores the outputted data. Device 1804 may also be a communication device to a manufacturing machine such as a CNC machine or any other CAM machine.
Device 1804 may also be a display that displays the generated CAM data in relation to a medical image. The processor also has a communication channel 1807 to receive external data from a communication device and to transmit data to an external device. The system in one embodiment of the present invention has one or more input devices 1805, which may be a keyboard, a mouse or any other device that can generated data to be provided to processor 1803. The processor can be dedicated hardware. However, the processor can also be a CPU or any other computing device that can execute the instructions of 1802. Accordingly, the system as illustrated in FIG. 32 provides a system for data processing resulting from a sensor such as medical image data or any other data source and is enabled to execute the steps of the methods as provided herein as an aspect of the present invention.
The methods as provided herein are, in one embodiment of the present invention, implemented on a system or a computer device. FIG. 33 illustrates a system 2600 used in generating CAM model of an implant in accordance with an aspect of the present invention. A processor 2604 is programmed to perform the steps of a method provided herein in accordance with an aspect of the present invention. On an input 2602 data related to a medical image of a patient, such as a patient bone such as a distal femur, is received. The processor processes the image data as described above to provide CAM model of a medical implant with a match between the generated medical implant and the patient's bone and outputs a data file representing the CAM model of the medical implant on 2605. Output 2605 may be connected to a manufacturing system 2606 that is enabled to provide the actual medical implant in its appropriate shape. For instance the CAM model may be provided or transformed into CNC instructions to a CNC machine that mills a material into the medical implant. In a next step the medical implant generated from the CAM model is removed from the manufacturing machine 2606 and implanted in a patient 2610 for instance by a surgeon.
There may be a choice of material for the implant that can be selected or the selected implant will be manufactured in one pre-determined material. In any event, in one embodiment of the present invention data of an implant related to patient image data is provided with manufacturing data related to manufacturing system 2606. For instance, if 2606 includes a CNC machine, then data related to cutting tools, preferred milling speeds and other data related to machining a piece of a pre-determined material is attached to the CAM file, so that the set-up of the manufacturing machine can be done based on data that is part of the customized CAM file. In a further embodiment of the present invention, an order of milling steps by manufacturing machine 2606 may be done based on an analysis of the customized CAM file. In a further embodiment an order of at least two milling or cutting steps may be pre-set in the CAM file, including a switching of tools. Availability of such data will minimize the need for pre-manufacturing planning and human intervention during manufacturing.
In one embodiment of the present invention an image storage device or a medical imaging device, the system 2604 and the manufacturing system 2606 are all connected via a network. In a further embodiment the network is the Internet. In yet a further embodiment of the present invention the system 2604 is authorized to receive data 2602 and 2606 is authorized to receive data 2605 from 2604 over the Internet.
In one embodiment of the present invention, the implant is already manufactured in oversized format, and wherein customization includes the removal of excess material to achieve the final and customized dimensions of the custom implant. This allows a core of an implant, which in many cases is the same for different patients, to be pre-manufactured. This allows for a fine-tuning of a custom implant that in many cases is faster than building an implant from scratch. The manufacturing machine to that purpose is provided with high accuracy measuring equipment that determines the actual dimensions of a pre-formed implant 2609 as shown in FIG. 33. Based on the received CAD or CAM file the pre-formed implant 2609 is milled to the exact measurements required for the custom implant.
In one embodiment of the present invention a manufacturing machine is a deposit machine which deposits material based on instructions which are derived from the CAD or CAM file. In that case a pre-form is a skeleton structure that is under-sized compared to the final dimensions of the custom implant. The manufacturing machine is equipped with high accuracy measuring devices. Material is deposited onto the skeleton or preform and measured until the desired dimensions have been achieved.
The finalized custom implant is removed from the manufacturing machine. It may undergo additional treatment, including finalizing treatment such as annealing, hardening, polishing, sterilizing testing, marking or any other treatment that is required to prepare for surgical insertion. The custom implant is provided to a surgeon or a surgical robot in an operating room and is implanted in the patient.
Thus, methods and systems are provided to create a customized medical implant or a patient specific medical implant based on a medical image of the patient and to create a CAD model of an implant that corresponds to an image of a bone or an organ of the patient, and to provide the CAD file to a manufacturing machine to manufacture the patient specific medical implant.
In order to further illustrate some of the steps of the present invention, triangulated figures have been enlarged to show a more detailed figure. Items 2601 and 2602 are shown enlarged in FIGS. 34-37. FIG. 34 illustrates an enlarged item 2601 which now shows the triangulation with a part 3401 which is shown in further enlargement in FIG. 35. FIG. 36 shows an enlarged 2602 and FIG. 37 shows an enlarged part 3601 in 2602. FIG. 38 shows element 2701 of which an enlarged part is shown in FIG. 39. FIG. 40 shows enlarged 2702 from FIG. 27 and part 4001 is further enlarged in FIG. 41. FIG. 42 shows part 7703 from FIG. 27 and part 4201 is further enlarged in FIG. 43.
In summary, in accordance with one aspect of the present invention, a method of making an implant from an image of a bone is provided. A processor receives image data of the image of the bone, segments a desired bone from the image data to create a segmentation of the desired bone, creates a polygonal mesh from the segmentation of the desired bone, extracts a part of the polygonal mesh corresponding to an implant region and parameterizes the part of the polygonal mesh so that every vertex therein is associated with two parameter values. The processor also resamples the polygonal mesh to generate a grid of 3D points regularly spaced in parametric domain and fits a B-spline surface to the grid using a least squares process, which can be weights. The processor can also smooth the surface.
Also in summary, in accordance with another aspect of the present invention, a method of making a patient specific implant is provided. In accordance with this method, an image of a portion of a patient's bone is taken. A processor is used to select a portion of the image to form a selected part of the image and processes the selected part of the image to form a surface. The processor creates a CAD model from the surface. Then a CAM system is used to manufacture the patient specific implant based on the CAD model. The method also comprising implanting the patient specific implant in the patient. In accordance with another aspect of the present invention, the CAM system can manufacture the customized medical implant from a preformed implant.
A solid CAD model can be created from the B-spline surface. Assembly features can be added to the CAD model. The assembly features can be a pin and screw thread. The implant can be manufactured using a CAM system. The polygonal mesh can be created using a Marching Cubes process or an Afront process or any other applicable process.
The following references provide background information generally related to the present invention and are hereby incorporated by reference: [1] W. Lorensen and H. Cline. “Marching cubes: A high resolution 3d surface construction algorithm,” in Proceedings of the 14th annual conference on Computer graphics and interactive techniques. ACM. 1987. p. 169; [2] J. Schreiner. C. Scheidegger. and C. Silva. “High-quality extraction of isosurfaces from regular and irregular grids,” IEEE Transactions on Visualization and Computer Graphics, pp. 1205-1212, 2006; [3] M. S. Floater and K. Hormann, “Surface parameterization: a tutorial and survey,” in In Advances in Multiresolution for Geometric Modelling. Springer, 2005, pp. 157-186; [4] A. Sheffer, E. Praun, and K. Rose, “Mesh parameterization methods and their applications,” in Foundations and Trends in Computer Graphics and Vision. Now Publishers, 2006, pp. 105-171; [5] M. O. Carmo, Differential Geometry of Curves and Surfaces. Prentice Hall, 1976; [6] E. Cohen, R. Riesenfeld, and G. Elber, Geometric modeling with splines: an introduction. AK Peters Ltd, 2001; [7] L. Piegl and W. Tiller, The NURBS book. Springer Verlag, 1997; [8] C. De Boor, A practical guide to splines. Springer Verlag, 2001; [9] M. Botsch, M. Pauly, L. Kobbelt, P. Alliez, B. Levy, S. Bischoff, and C. Rossi, “Geometric modeling based on polygonal meshes,” in SIGGRAPH '07: ACM SIGGRAPH 2007 courses. New York, N.Y., USA: ACM, 2007; [10] U. Pinkall and K. Polthier, “Computing discrete minimal surfaces and their conjugates,” Experimental Mathematics, vol. 2, pp. 15-36, 1993; [11] V. Weiss, L. Andor, G. Renner, and T. Varady, “Advanced surface fitting techniques,” Computer Aided Geometric Design, vol. 19, no. I, pp. 19-42, 2002; and [12] S. Open Cascade, “Open cascade technology ver 6.3,” 2010, URLwww.opencascade.org.
While there have been shown, described and pointed out fundamental novel features of the invention as applied to preferred embodiments thereof, it will be understood that various omissions and substitutions and changes in the form and details of the methods and systems illustrated and in its operation may be made by those skilled in the art without departing from the spirit of the invention. It is the intention, therefore, to be limited only as indicated by the scope of the claims.